﻿#include "noah.h"
#define PI 3.141592653589793238462643383f
#define WALL_BUFFER	0.00001f					//for collision detection
#define K 1.0
#define THRESHOLD 0.8f

/**
Temperature: 295.15K
Density of water at 295.15K: 0.9977735 g/mL


* K is dependent on the volume of the simulation.
Since pV = nRT where n is the number of particles in mol we can consider the right side
constant, K. R is the ideal gas constant and T is temperature. So temperature could be modeled
if we wanted to get fancy, but for now we assume room temperature of 295.15Kelvin and R = 8.314472 m^2 kg^s-2 Kelvin-1
so after finding n we can solve for K.

n = Volume(L) * density(kg/L) / molecular mass(kg/mol)

****
note we need to calculate Volume if we want to add particles
****
Molucular mass of Water = 18.01528 g/mol = 0.01801528 kg/mol

Assuming 10mL (Meaning for 100 particles each particle is 0.1mL), ***** Need to adjust the mass of the particles appropriately*****

n = 0.01(L) * 1(kg/L) / 0.01801528 kg/mol
n = 0.55508435061791989910786843168688

so,

K = nRT
= 0.55508435061791989910786843168688(mol) * 295.15(K) * 8.314472 (m^2/K*kg^s2)
K = 1362.1861057946365529705894107669

*/


struct WindowSize {
	int widthPixel;
	float width;
	int heightPixel;
	float height;
	float zNear;
	float zFar;

	WindowSize() :
	widthPixel(512), heightPixel(512),
		width(20.0f), height(20.0f),
		zNear(1.0f), zFar(20.0f)
	{}
};

int g_particleIndex = 0;
WindowSize g_winSize;

// Integration
enum IntegrationMethods {
	EULERIAN,
	LEAPFROG
};
IntegrationMethods currentIntegrationMethod = LEAPFROG;

// GLUI
int   main_window;
GLUI* glui;

bool leftMouseButtonDown = false;
bool rightMouseButtonDown = false;
vec2f g_mouseCoord;

//Precomputed values
float g_gradient_spiky;
float g_6poly;
float g_gradient_6poly;
float g_laplacian_6poly;
float g_laplacian_viscosity;

//Particle jazz
std::vector<Particle> g_particles;
std::map<int, GLfloat> g_particleDensity;
std::map<int, vec3f> g_pressureForce;
std::map<int, vec3f> g_viscosityForce;
std::map<int, vec3f> g_colourFieldGradient;
std::map<int, GLfloat> g_colourFieldLaplacian;

//Simulation Globals
GLfloat g_particleMass; // kg
GLfloat g_smoothingLength = 0.01f;
GLfloat g_particleVolume; // m^3
GLfloat g_restDensity;  // kg/m^3
GLfloat g_viscosityCoeff; 	// pascal-second (Pa.s) = 1 kg m^-1 s^-1
GLfloat g_surfaceTensionCoeff;
GLfloat g_simScale = 0.004f;
GLfloat g_particleRadius;
vec3f g_gravity(0.0f, -0.1f, 0.0f);
GLfloat g_timeStep;
int g_enable3d = 1;
GridMode gridMode = HIDDEN;
bool filledGrid = false;

// For faster rendering
GLuint sphereDisplayList;

// For collision detection
GLfloat g_stiff = 10000.0;
GLfloat g_damp = 256.0;

GLfloat g_initParticleSpacing;
GLfloat g_scaledSpacing;
vec3f g_volMin(-10.0f, -10.0f, -10.0f);
vec3f g_volMax(10.0f, 10.0f, 10.0f);
Grid g_grid(g_volMin, g_volMax, 0.007f / g_simScale);

vec3f g_initParticlePosMin(-9.0f, -7.0f, -7.0f);
vec3f g_initParticlePosMax(-4.0f, 9.0f, 9.0f);

vec3f g_volDimensions( g_volMax.x() - g_volMin.x(), g_volMax.y() - g_volMin.y(), g_volMax.z() - g_volMin.z() );

//Camera variables
GLfloat g_camX = 0.0f;
GLfloat g_camY = 0.0f;
int g_prevX = 0;
int g_prevY = 0;

//Values to be precomputed
GLfloat pre_w_poly6; 
GLfloat smoothing_len_squared;
GLfloat pre_lapacian_w_poly6;
GLfloat pre_gradient_w_poly6;
GLfloat pre_gradient_w_spiky;
GLfloat pre_laplacian_w_viscosity;

void reshape(int width, int height) 
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	GLfloat minDim = std::min(g_winSize.width, g_winSize.height);
	if(width > height) 
	{
		g_winSize.width = minDim * (GLfloat)width/(GLfloat)height;
		g_winSize.height = minDim;
	} else 
	{
		g_winSize.height = minDim * (GLfloat)height/(GLfloat)width;
		g_winSize.width = minDim;
	}

	g_winSize.widthPixel = width;
	g_winSize.heightPixel = height;
	// Allows the water container to resize, while preserving the size of the particles.
	g_volMax.d_comp[0] = width / 50;
	g_volMin.d_comp[0] = -width / 50;
	g_volDimensions = vec3f( g_volMax.x() - g_volMin.x(), g_volMax.y() - g_volMin.y(), g_volMax.z() - g_volMin.z() );
	glViewport( 0, 0, g_winSize.widthPixel , g_winSize.heightPixel );
	g_grid.resize(g_volMin, g_volMax);
	gluPerspective( 60.0f, g_winSize.width / (g_winSize.height), 10.0f, 800.0f);
}

void mouse(int button, int state, int x, int y) 
{
	if (button == GLUT_LEFT_BUTTON)
	{
		leftMouseButtonDown = !leftMouseButtonDown;
	} else if (button == GLUT_RIGHT_BUTTON) {
		rightMouseButtonDown = !rightMouseButtonDown;
	}
	g_prevX = x;
	g_prevY = y;
}

void mouseMotion(int x, int y) 
{
	g_mouseCoord.d_comp[0] = static_cast<float>(x);
	g_mouseCoord.d_comp[1] = static_cast<float>(y);
	if (rightMouseButtonDown) {
		int deltaX = x - g_prevX;
		int deltaY = y - g_prevY;

		g_camX += deltaX;
		g_camY += deltaY;

		g_prevX = x;
		g_prevY = y;
	}
}

void precompute()
{
	pre_w_poly6 = (315.0f/(64.0f*PI*(std::pow(g_smoothingLength,9)))); 
	smoothing_len_squared = g_smoothingLength*g_smoothingLength;
	pre_lapacian_w_poly6 = (945.0f / (8.0f*PI*std::pow(g_smoothingLength, 9)));
	pre_gradient_w_poly6 = -1 * (945.0f/(32.0f*PI*std::pow(g_smoothingLength, 9)));
	pre_gradient_w_spiky = -1 * (45.0f/(PI*std::pow(g_smoothingLength, 6)));
	pre_laplacian_w_viscosity = (45.0f/(PI*std::pow(g_smoothingLength, 6)));
}

GLfloat W_poly6(float length, GLfloat h)
{
	return pre_w_poly6 * std::pow((smoothing_len_squared - length*length), 3);
}

GLfloat laplacian_W_poly6(GLfloat length, GLfloat h)
{
	GLfloat length_squared = length*length;
	return pre_lapacian_w_poly6 * (smoothing_len_squared - length_squared) * (length_squared - (0.75f * (smoothing_len_squared - length_squared)));
}

vec3f gradient_W_poly6(vec3f r, GLfloat h)
{
	float length = r.length(); 
	r *= pre_gradient_w_poly6  * std::pow((smoothing_len_squared - length*length), 2);
	return r;
}

vec3f gradient_W_spiky(vec3f r, GLfloat h)
{
	float length = r.length(); 
	r *= pre_gradient_w_spiky/length*std::pow((g_smoothingLength - length), 2);
	return r;
}

GLfloat laplacian_W_viscosity(GLfloat length, GLfloat h)
{
	return pre_laplacian_w_viscosity*(g_smoothingLength-length);
}

void update(){
	if (leftMouseButtonDown) 
	{
		GLfloat xCoord = (g_mouseCoord.d_comp[0] / (GLfloat)g_winSize.widthPixel) * (g_volMax-g_volMin).x() + g_volMin.x();
		GLfloat yCoord = (g_mouseCoord.d_comp[1] / (GLfloat)g_winSize.heightPixel) * (g_volMax-g_volMin).y() + g_volMin.y();
		g_particles.push_back(Particle(vec3f(xCoord,-yCoord, 0.0f), vec3f(1.0f, -0.4f, 0.03f)));
	}
	g_grid.update(g_particles);
	GLfloat h = g_smoothingLength;

	g_pressureForce.clear();
	g_viscosityForce.clear();
	g_colourFieldGradient.clear();
	g_colourFieldLaplacian.clear();
	g_particleDensity.clear();

	// calculate densities
	for (std::vector<Particle>::iterator particle = g_particles.begin(); particle != g_particles.end(); ++particle)
	{
		std::vector<Particle*> neighbours = g_grid.findNeighbours(particle->pos);
		for (std::vector<Particle*>::iterator it = neighbours.begin(); it != neighbours.end(); ++it)
		{
			Particle* nParticle = *it;
			if (nParticle->id != particle->id) {
				vec3f r = (particle->pos - nParticle->pos) * g_simScale;
				GLfloat length = r.length();
				if(length <= h)
				{
					g_particleDensity[particle->id] += g_particleMass * W_poly6(length, h);
				}
			}
		}
	}

	// calculate forces and colour-field
	for (std::vector<Particle>::iterator particle = g_particles.begin(); particle != g_particles.end(); ++particle)
	{
		std::vector<Particle*> neighbours = g_grid.findNeighbours(particle->pos);
		for (std::vector<Particle*>::iterator it = neighbours.begin(); it != neighbours.end(); ++it)
		{
			Particle* nParticle = *it;
			if (nParticle->id != particle->id) {
				vec3f r = (particle->pos - nParticle->pos) * g_simScale;

				GLfloat length = r.length();
				if(length <= h)
				{
					GLfloat particleDensity = g_particleDensity[particle->id];
					GLfloat neighbourDensity = g_particleDensity[nParticle->id];

					GLfloat particlePressure = 0.5f * (particleDensity - g_restDensity);
					GLfloat neighbourPressure = 0.5f * (neighbourDensity - g_restDensity);

					//Symmetric Pressure, adheres to Newton's 3rd law
					g_pressureForce[particle->id] -= gradient_W_spiky(r, h) * 0.5f * g_particleMass * ((particlePressure + neighbourPressure) / neighbourDensity);
					g_viscosityForce[particle->id] += (nParticle->vel - particle->vel) * g_viscosityCoeff * (g_particleMass / neighbourDensity) * laplacian_W_viscosity(length, h);
					g_colourFieldGradient[particle->id] += gradient_W_poly6(r, h) * (g_particleMass / neighbourDensity);
					g_colourFieldLaplacian[particle->id] +=  (laplacian_W_poly6(length, h) * g_particleMass /  neighbourDensity);
				}
			} 
		}
	}

	float speedLimit = 200;
	float speedLimitSquared = speedLimit * speedLimit;
	// move particles
	for (std::vector<Particle>::iterator particle = g_particles.begin(); particle != g_particles.end(); ++particle)
	{
		GLfloat gradientLength = g_colourFieldGradient[particle->id].length();
		vec3f surfaceTensionForce(0.0f, 0.0f, 0.0f);
		if(gradientLength >= THRESHOLD) // (2.29)
		{ 	
			surfaceTensionForce = (g_colourFieldGradient[particle->id] *(1 / gradientLength)) * g_colourFieldLaplacian[particle->id] * (-1 * g_surfaceTensionCoeff);
		}

		vec3f totalForce = g_pressureForce[particle->id] + g_viscosityForce[particle->id] + surfaceTensionForce;

		vec3f acceleration;
		if (g_particleDensity[particle->id]) {
			acceleration = (totalForce  * (1/g_particleDensity[particle->id]));
		}

		
		/*
		 Collision Detection & Leap Frog integration based off of Fluid v.2 simulation found 
		 http://www.rchoetzlein.com/eng/graphics/fluids.htm
		*/
		// Velocity limiting 
		float speed = acceleration.dot(acceleration);

		if ( speed > speedLimitSquared ) {
			acceleration *= speedLimit / sqrt(speed);
		}

		// X-axis walls
		float diff = 2 * g_particleRadius - (particle->pos.x() - g_volMin.x()) * g_simScale;
		float adj;
		if (diff > WALL_BUFFER) {
			adj = g_stiff * diff - g_damp * particle->vel.x();
			acceleration.x() += adj;
		}
		diff = 2 * g_particleRadius - (g_volMax.x() - particle->pos.x()) * g_simScale;
		if (diff > WALL_BUFFER) {
			adj = g_stiff * diff - g_damp * -particle->vel.x();
			acceleration.x() += -adj;
		}

		//Y-axis walls
		diff = 2 * g_particleRadius - (particle->pos.y() - g_volMin.y()) * g_simScale;
		if (diff > WALL_BUFFER) {
			adj = g_stiff * diff - g_damp * particle->vel.y();
			acceleration.y() += adj;
		}

		diff = 2 * g_particleRadius - (g_volMax.y() - particle->pos.y()) * g_simScale;
		if (diff > WALL_BUFFER) {
			adj = g_stiff * diff - g_damp * -particle->vel.y();
			acceleration.y() += -adj;
		}

		// Z-axis walls
		diff = 2 * g_particleRadius - (particle->pos.z() - g_volMin.z()) * g_simScale;
		if (diff > WALL_BUFFER) {
			adj = g_stiff * diff - g_damp * particle->vel.z();
			acceleration.z() += adj;
		}
		diff = 2 * g_particleRadius - (g_volMax.z() - particle->pos.z()) * g_simScale;
		if (diff > WALL_BUFFER) {
			adj = g_stiff * diff - g_damp * -particle->vel.z();
			acceleration.z() += -adj;
		}

		// Add external forces
		acceleration += g_gravity;

		switch (currentIntegrationMethod) {
			case LEAPFROG:
			{
				vec3f vnext = acceleration;							
				vnext *= g_timeStep;
				vnext += particle->vel;						// v(t+1/2) = v(t-1/2) + a(t) dt
				particle->vel = particle->vel;
				particle->vel += vnext;
				particle->vel *= 0.5;					// v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5		used to compute forces later
				particle->prevVel = vnext;
				vnext *= g_timeStep/ g_simScale;
				particle->pos += vnext;						// p(t+1) = p(t) + v(t+1/2) dt
				
				break;
			}
			case EULERIAN:
			{
				
				particle->vel += acceleration * g_timeStep;
				particle->pos += particle->vel * g_timeStep;
				break;
			}
			default:
				break;
		}
		if (!g_enable3d) {
			particle->pos.d_comp[2] = 0.0f;
		}
	}
}

int frame = 0;
float timebase, currenttime = 0;
int fps = 0;
char titleFPS[20];
void render() 
{
	frame++;
	if (frame % 100 == 0) {
	currenttime = static_cast<float>(glutGet(GLUT_ELAPSED_TIME));
		fps = ceil(frame*1000.0/(currenttime-timebase));
		sprintf(titleFPS, "%i fps", fps);
		glutSetWindowTitle(titleFPS);
		timebase = currenttime;		
		frame = 0;
	}

	update();
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glColor3f(0.5f,0.5f,0.5f);

	glTranslatef(0.0f, 0.0f, -35.0f);
	glRotatef(g_camY,1.0f,0.0f,0.0f);
	glRotatef(g_camX,0.0f,1.0f,0.0f);

	glPushMatrix();
	glScalef(g_volDimensions.x(), g_volDimensions.y(), g_volDimensions.z());
	glutWireCube(1.0f);
	glPopMatrix();

	glEnable(GL_LIGHTING);
	glEnable ( GL_LIGHT0 );
	glEnable(GL_COLOR_MATERIAL);

	GLfloat diff[4];
	GLfloat spec[4];
	GLfloat shininess = 10.0;

	diff[0] = 0.0f; diff[1] = 0.0f; diff[2] = 0.8f; diff[3] = 1.0f;
	spec[0] = 0.1f; spec[1] = 0.1f; spec[2] = 0.1f; spec[3] = 0.1f;
	glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, &diff[0]);

	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, &spec[0]);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, &shininess);		
	glColorMaterial( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE );
	// If the user is displaying the grid as filled in cubes, then don't show particles.
	if (gridMode != MINECRAFT) {
		for (std::vector<Particle>::iterator it = g_particles.begin(); it != g_particles.end(); ++it) {
			glPushMatrix();
			glColor4f(0.0f, 0.5f, 1.0f, 0.5f);
			glTranslatef(it->pos.x(), it->pos.y(), it->pos.z());
			glCallList(sphereDisplayList);
			glPopMatrix();
		}
	}
	//Time to render the digital frontier
	if(gridMode != HIDDEN)
	{
		g_grid.drawGrid(gridMode);
	}
	glDisable(GL_LIGHTING);

	glFlush();
	glutSwapBuffers();
	glutPostRedisplay();
}

void myGlutIdle( void )
{
}

void initParticles() {
	g_particleVolume = 0.00000030543f;
	g_restDensity = 600.0f;
	g_particleMass = g_particleVolume * g_restDensity;
	g_viscosityCoeff = 0.2f;	// pascal-second (Pa.s) = 1 kg m^-1 s^-1
	g_surfaceTensionCoeff = 0.0000728f;
	g_gravity = vec3f(0.0f, -9.8f, 0.0f);
	g_timeStep = 0.003f;
	g_particleRadius = 0.004f;
	g_initParticleSpacing = powf(g_particleMass / g_restDensity, 1.0f/3.0f);
}

void createParticles( ) {
	GLfloat scaledSpacing = g_initParticleSpacing / g_simScale;	

	for (float y = g_initParticlePosMax.y(); y >= g_initParticlePosMin.y(); y-= scaledSpacing) {
		for (float z = g_initParticlePosMin.z(); z <= g_initParticlePosMax.z(); z+= scaledSpacing) {
			for (float x = g_initParticlePosMin.x(); x <= g_initParticlePosMax.x(); x+= scaledSpacing) {
				g_particles.push_back(Particle(vec3f(x,y,z), vec3f(0, 0, 0)));
			}
		}
	}
}

void reset()
{
	g_enable3d = 1;
	initParticles();
	g_particles.clear();
}

void buttonHandler(int control) 
{
	if (control == 0) 
	{
		createParticles();
	} 
	else if (control ==1) 
	{	
		reset();
	}
	glui->sync_live();
}

void keyboard(unsigned char c, int x, int y) 
{
	//Dam Break
	if (c == '1')
	{
		reset();
		initParticles();

		g_initParticlePosMin = vec3f(-9.0f, -7.0f, -7.0f);
		g_initParticlePosMax = vec3f(-4.0f, 9.0f, 9.0f);

		createParticles();
	}
	else if (c == '2') //Egyptian Killer
	{
		reset();
		initParticles();

		g_initParticlePosMin = vec3f(-9.0f, -10.0f, -7.0f);
		g_initParticlePosMax = vec3f(-4.0f, 10.0f, 9.0f);

		createParticles();

		g_initParticlePosMin = vec3f(4.0f, -10.0f, -7.0f);
		g_initParticlePosMax = vec3f(9.0f, 10.0f, 9.0f);

		createParticles();
	}
	else if (c == '3') //Hyper Cube
	{
		reset();
		initParticles();
		g_initParticleSpacing /= 1.25f;
		g_initParticlePosMin = vec3f(-5.0f, 0.0f, -5.0f);
		g_initParticlePosMax = vec3f(5.0f, 10.0f, 5.0f);

		createParticles();
	}
	else if (c == '4') //Power Tower
	{
		reset();
		initParticles();
		g_viscosityCoeff = 3.5f;
		g_initParticlePosMin = vec3f(-3.5f, 10.0f, -3.5f);
		g_initParticlePosMax = vec3f(3.5f, 50.0f, 3.5f);

		createParticles();
	}
	else if (c == '5') //PopCorn
	{
		reset();
		g_viscosityCoeff = 0.7021f;
		g_surfaceTensionCoeff = 0.2515f;
		g_initParticlePosMin = vec3f(-3.0f, 0.0f, -3.0f);
		g_initParticlePosMax = vec3f(3.0f, 10.0f, 3.0f);

		createParticles();
	}
	else if (c == 'e') 
	{
		currentIntegrationMethod = EULERIAN;
		g_timeStep = 0.008f;
	}
	else if (c == 'l') 
	{
		currentIntegrationMethod = LEAPFROG;
		g_timeStep = 0.003f;
	}
	else if(c == 'g')
	{
		gridMode = (GridMode)((gridMode + 1) % 4);
	}
}

void init() 
{
	glutDisplayFunc(render);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard);
	glutMouseFunc(mouse);
	glutIdleFunc(myGlutIdle);
	glutMotionFunc(mouseMotion);

	glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
	glEnable(GL_DEPTH_TEST);

	glMatrixMode(GL_PROJECTION);

	glLoadIdentity();

	glFrustum( -g_winSize.width/2.0, g_winSize.width/2.0,
		-g_winSize.height/2.0, g_winSize.height/2.0,
		g_winSize.zNear, g_winSize.zFar);

	// Create the display list
	sphereDisplayList = glGenLists(1);
	glNewList(sphereDisplayList, GL_COMPILE);
	glutSolidSphere(0.6, 20, 20);
	glEndList();
	
	std::cout<<"[1-5] shifts between initial starting states"<<std::endl;
	std::cout<<"E - Eulerian integration"<<std::endl;
	std::cout<<"L - Leapfrog integration(default"<<std::endl;
	std::cout<<"g - Cycle through grid modes"<<std::endl;
}

int main(int argc, char *argv[]) 
{  
	// let glut eat any command line args that it owns
	glutInit( &argc, argv );
	// full color, double buffered
	glutInitDisplayMode( GLUT_RGBA | GLUT_DEPTH | GLUT_DOUBLE );
	glutInitWindowSize(g_winSize.widthPixel, g_winSize.heightPixel);
	glutCreateWindow( "NOAH" );
	glui = GLUI_Master.create_glui( "Settings" );
   
	new GLUI_Button(glui, "Reset", 1, buttonHandler);
	new GLUI_Button(glui, "Add Water Cube", 0, buttonHandler); 
	(new GLUI_Spinner( glui, "Rest Density:", &g_restDensity ))
		->set_float_limits( 0.00001f, 1000.0f );
	(new GLUI_Spinner( glui, "Viscosity Coeff:", &g_viscosityCoeff ))
		->set_float_limits( 0.00001f, 10.0f );
	(new GLUI_Spinner( glui, "Surface Tension:", &g_surfaceTensionCoeff ))
		->set_float_limits( 0.00001f, 10.0f );
	(new GLUI_Spinner( glui, "Gravity:", &g_gravity.y() ))
		->set_float_limits( -10.0f, 10.0f );
	(new GLUI_Spinner( glui, "Timestep:", &g_timeStep ))
		->set_float_limits( 0.0003f, 0.003f );
	new GLUI_Checkbox(glui, "3d mode:",&g_enable3d,2,buttonHandler);

	glui->set_main_gfx_window( main_window );

	/* We register the idle callback with GLUI, *not* with GLUT */
	GLUI_Master.set_glutIdleFunc( myGlutIdle ); 

	// setup remaining state
	init();

	precompute();

	initParticles();
	createParticles();
	glui->sync_live();

	// start processing events
	glutMainLoop();

	return 0;
}


